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O ■ Abstract 

In order to quantitatively study the accuracy of the unconditionally stable coars- 
ening algorithms, we calculate the Fourier space multi step error on the order pa- 
rameter field by explicitly distinguishing the analytic time r and the algorithmic 
time t. The calculation determines the error in the order parameter and the scaled 
•th ■ correlations. This error contributes a correction term in the analytic time step, 

£>y which is crucial in understanding the accuracy in unconditionally stable coarsening 

i-Q , algorithms. 
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O ■ 1 Introduction 
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The Partial Differential Equations (PDEs) that governs conserved and non- 
conserved scalar coarsening systems are the Cahn-Hilliard (CH) p] and Allen- 
Cahn (AC) equations [2] respectively. Examples are found in polymer mixtures 
[3], alloys [4T5] . liquid- crystals [5117] . and in cosmology [S]. These equations, as 
well as their various extended forms, are an active area of research. They model 
the phase separation that occurs after a quench from a high temperature dis- 
ordered phase into two distinct phases at low temperatures. The pattern of the 
two phase regions coarsens as the time increases, i.e., the length-scale of these 
regions grows. The scaled correlation is an important probe to characterize 
the late time scaling regime, where the dynamics are dominated by a single 
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length scale, the average domain size L, which increases with a power law in 
time r, L(t) ~ r a [9j. For the scalar order parameter systems considered in 
this paper, a = 1/3 and 1/2 for conserved and and non-conserved systems, re- 
spectively 0. To study the coarsening systems, one needs to choose a random 
initial condition, which corresponds to a high temperature disordered phase. 
Since it is essentially not possible to find an exact solution for random ini- 
tial conditions, computer simulation plays an important role in understanding 
these systems. 

Because of the nature of slow coarsening process in the late time scaling regime, 
the main challenge in computer simulation of coarsening systems is to develop 
an efficient algorithm. The most straightforward approach is the Euler algo- 
rithm, which must employ a time step At^ u ~ (Ax) 4 and (Ax) 2 to maintain 
stability for conserved and non-conserved systems, respectively, where Ax is 
the lattice spacing. On the other hand, one has to use a lattice spacing Ax < £ 
to resolve the interfacial profile, where £ is the interfacial width. This implies 
that the Euler update is inefficient, and in practice computationally costly to 
use to evolve large systems until the late time scaling regime. Various com- 
putational algorithms have been developed by increasing the algorithmic time 
step At compared to the simplest Euler discretization. For example, the Cell 
Dynamical System (CDS) pT0]llllll2j was broadly applied [TBim] as an efficient 
technique. However, a transformation of CDS to PDE pT5]fT6] showed that CDS 
changes the form of the equation of motion within the same universality class 
and gets a less strict stability threshold for its Euler update, but the Euler 
time step is still limited to (Ax) 4 and (Ax) 2 for conserved and non-conserved 
systems, respectively. More recently, conditionally stable algorithms such as 
Fourier spectral methods [TTfT8] have been developed by introducing implicit 
terms in the update equation, and have been shown to allow an increase of 
the algorithmic time step At by two orders of magnitude. However, none of 
these methods can adjust to the naturally slowing dynamics - they eventually 
become more and more inefficient at the late time scaling regime. 

The recently developed unconditionally stable algorithms |19|l20f 21] overcomes 
the difficulty of a fixed time step for maintaining stability and allows for 
arbitrarily large algorithmic time steps. In CH equation, the algorithm allows 
an optimal driving scheme At = At 2 / 3 with the error in correlations scales 
as yA for small prefactor A [22|23] . This implies that any driving scheme 
that is slower than the optimal one is wastefully accurate because the error is 
negligible in correlations [23]. This algorithm has also been applied to other 
interesting systems, such as the Phase Field Crystal (PFC) model [2Hl2"o] 
and the Swift-Hohenberg (SH) equation and showed significantly improved 
efficiency |26j . 

While we have no doubt about the efficiency of unconditionally stable algo- 
rithms, the accuracy of unconditionally stable algorithm, remains somewhat 



mysterious, and some open questions remain unanswered. For example, why 
does the single step error accumulate over time [22|23] ? why does the fitting 
for the maximal analytic time step in non-conserved system get worse as a [see 
Eq. (12~4|) for a precise definition] increases? In this paper, by explicitly calcu- 
lating the Fourier space multi step error using Eq. (fT5l) . we are able to answer 
these questions. The rest of the paper is organized as follows. In Section 2, we 
briefly review the results in the previous work [19. 20. 21. 22.23] and introduces 
the basic concepts and notations of unconditionally stable algorithms. In Sec- 
tion 3, we study the Fourier space multi step error, where the distinction of 
the analytic time and the algorithmic time is necessary. Based on this result, 
in Section 4, we determine the accuracy in the scaled structure factor g(x) 
and h(x). We then calculate the analytic time step in Section 5, taking into 
account the error in h(x). The summary are presented in Section 6. In the 
appendix, we derive the scaling of field derivatives for all modes in Fourier 
space. For simplicity but without loss of generality, we restrict our analysis to 
two dimensions (2D). 



2 Review of Unconditionally stable algorithms 



For coarsening systems, a suitable free energy functional to describe the or- 
dered phase is [9] 



d x 



i|V0| 2 + V(0) 



where 0(x, t) is the order parameter, and the potential V{4>) has a double- well 
structure V{4>) = (0 2 — l) 2 /4. The two minima of V correspond to the two 
equilibrium phases <fi = ±1. The term V(4>) describes the energy in the bulk. 
The gradient- squared term in Eq. ([1]) associates an energy cost to an interface 
between the two phases. 

When the order parameter is conserved under the dynamics, the equation of 
motion can be written in the form of a continuity equation, d<fi/dt = — V • 
j, with current j = — MV(<5.F/<50), where M is the mobility. Absorbing M 
into the time scale, we obtain the dimensionless form of Cahn-Hilliard (CH) 
equation pQ: 



^ = V 2 ^ = -V 2 (V 2 + 0-0 3 ). (2) 

In the case where the order parameter is not conserved by the dynamics, the 
simplest dissipative equation for the time evolution of the field <p is 



Eq. (|3D describes that the rate of change of is proportional to the gradient 
of the free-energy functional in functional space, and is often called the Allen- 
Cahn (AC) equation [2]. 

In order to obtain unconditionally stable algorithms for these equations, we 
introduce a mix of explicit and implicit terms in the update [X91I20U21J . For 
the CH equation, we have 



fa+At - (oi - l)AtV 2 fa +A t - (o 2 - l)AtV 4 fa +A t 

= fa - At(a x V 2 fa + a 2 V 4 fa - V 2 t 3 ), (4) 

and for the AC equation, we have 

fa+At + (tti - l)Atfa +At + (02 - l)AtV 2 i+A i 

= t + At(ai0 4 + a 2 V 2 4 -0 i 3 ). (5) 

Implicit terms (fa+At) are on the left and the explicit terms (fa) are on the 
right. Unconditionally stable algorithms are obtained for a\ > 2 and a 2 < 1/2 
[T9|20|21j . In Fourier space, Eq. (J3J) can be directly solved as 

fa(t + At) = fa(t) + At eff (k,At)^-, (6) 

where the /c-dependent effective time step is 



Ateff{k > At) " l + Atk^a 1 -l) + (l-a 2 )^Y (7) 

and d<f)k/dr is the Fourier transform of d<f)/dr in Eq. (J2J) and is a function of 
<pk(t)- As Eq. (J2J) reveals, the unconditionally stable algorithm has a mode- 
dependent effective time step At e ff(k,Ai). The optimal driving scheme for 
CH equation is an algorithmic time step At = At 2 / 3 , where the prefactor A 
determines the accuracy [22|2"3] . 

On the other hand, the effective time-step for AC equation is 



At 
1 + At[(oi -Tf+ (1 - a2)k 2 



At eff (k,At) = .n-..,u,2\ 



which indicates that 



At eff (k, At) < At eff (k, oo) < -, (9) 

d\ — i 

where a\ — 1 > and 1 — 02 > are used. We cannot obtain an accelerated 
algorithm in non-conserved systems since there is an upper bound for the 
effective time-step At e ff(k, At). In fact, the optimal driving scheme for AC 
equation is an algorithmic time step At = 00 [22J. 



3 Fourier space multi step error 



In a computer simulation, assuming that we have updated the system m — 1 
times, the algorithmic time is t m _i, and the field is 4> k {t m _i) in Fourier space. 
We next update the system using an algorithmic time step At m , and the field 
evolves according to Eq. §§$): 

<Pk(t m ) = <pk{tm-l) + At e ff(k, At m )——, (10) 



where t m = £ m _i + A£ m . To obtain the error of this unconditionally stable algo- 
rithm A(pk, we need to compare the field evolved by an unconditionally stable 
algorithm 4> k {t) to the exact dynamics 0fc(r) evolved by the same amount of 
energy [22.23J. The analytic time step At and the analytic time r can be 
determined by the energy evolution of each update [23] . At an analytic time 
r m _i, the system evolution after an analytic time step Ar m is governed by the 
Taylor expansion: 

Mr m ) = Mr m -i) + E 7n?^f > (11) 



where r m = r m _i + Ar m . The Fourier space multi step error up to the mth 
update is 

A(j) k (m) = (j) k {t m ) - fc (r m ). (12) 

Note that 

A(j> k {m) = A(p k {m- 1) + 5(j) k (m), (13) 

where 



~ d n d) k AtI 



i*(m) = [At eff (k, At m ) - Ar m ] -^ - ^ -£± ~f ■ (14) 



We obtain that 



A(f) k (m) = ^<50 fc (m) « / 5<p k dm } (15) 

where we use the integral to approximate the summation. Eq. (fT5|) explains 
the numerical results that the single step error 8<\) k {m) accumulates over time, 
which was observed in previous studies [2~212"3"] . Note that as long as there is 
a distinction between the analytic time r (time step At) and the algorithmic 
time t (time step At), we can use Eq. (fT5l) to calculate the Fourier space multi 
step error. 

As an interesting application of the above result Eq. ffT51) . we study the ac- 
curacy in the simplest Euler algorithm. The Euler algorithm has a mode- 
independent fixed time step AtEu to update the system, and At e ff(k, At) = 
AtEu- In conserved systems, evolving from r to r max with a fixed time step 
AtEu ~ At ~ dr/dm [23J, we obtain the Fourier space multi step error as 
kL « 1: 



A0 fc ~- / / A ^V - 5 /^T 

V J At Eu 2 

TO 

_ 3At Eu / 2/3 -2/3N 3Atsu 2/3 

— 7 \^ T max T J ~ 7 T > l i0 J 

where we ignore the higher order terms in the sum and d 2 4>k/dr 2 ~ r -5 / 3 as 
fe ~ 1/L is used [23]. We can use this to determine the error in the scaled 
structure factor g(kL) = (<pk(p~k)/L 2 , where the angle brackets indicate an 
average over orientations and initial conditions. The absolute error is 



Ag(kL) 



2A0. 



kVk 



L 2 



^, (17) 



where <pk ~ L as k ~ 1/L is used [23]. Since L ~ r 1 / 3 — >• oo at late times, the 
error Ag(kL) — > at late times. This is why the Euler algorithm --or any 
fixed time step algorithm — is wastefully accurate. 



4 Error in scaled structures 



In this Section, we calculate the error in scaled structures due to the Fourier 
space multi step error determined in Section 3 with the optimal driving schemes 
discussed in Section 2. 



4-1 Conserved systems 



Using the formula discussed in Section 3, we obtain the Fourier space single 
step error with the optimal driving scheme At = Ar 2 ^ 3 [23] 



8h ~ ^^A^r-^k-^ (18) 

where A, B and C are constants [23], and dc^k/dr ~ r~ 5 / 6 A;~ 1//2 as k > 1/L is 
used [see Appendix]. For a small A, evolving from r to r max with time-step 
At = At 2 ! 3 w At ~ dr/dm, we obtain the Fourier space multi step error: 



^~ b ^ 1 ^ T-^ft-** 



LpCy/q! - \rj k -ili 1/6 (19) 

2R 1 / 3 'max- l ±J / 

We can use this result to determine the error in the scaled structure factor 
g{kL) = (0fc0_fc)/L 2 , where the angle brackets indicate an average over orien- 
tations and initial conditions. The absolute error is 



Ag(kL) 



2A0 



k<Pk 



L 2 



L Q Cy/ai - 1 \[A 

BV3 (A;L) 2 ' { ' 



where <pk ~ T mJx^ 3 ^ 2 as ^ ~ V-^ i s used [see Appendix]. We can subsequently 
obtain the absolute error in the scaled derivative structure factor h{kL) = 
(kL) 2 g(kL): 



Ah(kL) = (kL) 2 Ag(kL) ~ ^^^ VA. (21) 

To test Eq. (J2"U{) . we measure S(fc, t) = (4>k4>-k), and obtain the time-independent 
scaling function S'(y) = E 2 S{yE) using the energy density .E, where y = k/E. 
In Fig. [1] we plot the absolute difference between scaled structures obtained 
using an Euler algorithm and an unconditionally stable algorithm, AS(y) vs. 
y. We find that AS(y) ~ y~ 2 as y > 1, which implies Ag{x) ~ x~ 2 as x > 1, 
where x = fcL. We confirm the prediction in Eq. (1201) . 
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Fig. 1. The absolute difference between unconditionally stable and Euler updates, 
AS(y) vs. y = k/E for system size L^ = 512 at r = 1024. The Euler update used 
AtEu = 0.03 and the unconditionally stable update used At = O.Olr 2 ' 3 . The plot 
shows that AS(y) ~ y~~ 2 as y > 1, in agreement with Eq. (|20p . 

^.S Non-conserved systems 



Using the formula discussed in Section 3, we obtain the Fourier space single 
step error with the maximal driving scheme At = oo [22] 



(Xi — 1 



Ar„ 



r -3/4 A; -l/ 



-1/2 



(22) 



where d^^jdr ~ t _3 / 4 A; _1 / 2 as fe > 1/L is used [see Appendix]. Evolving 
from To to r max with time-step Ar ro ~ dr/dm, we obtain the Fourier space 
multi-step error: 



i ma. 



Ar n 



1 



T"0 

k- l ' 2 T l JU 



d\ — 1 

(7 



Ar rt 



3 / 4 ^ 



r ' or 



f tan" 1 (a/0 



(23) 



where 



a = 



1 - a 2 



V Oi — 1 
We then obtain the absolute error in scaled structure: 



(24) 



Ag(kL) 



L 2 



kVk 



(kLf 



£ tan- 1 (a/0 



-1 



(25) 



where <\> k ~ T raJx^ 3 ^ 2 as ^ ~ V-^ * s used [see Appendix]. We can subsequently 
obtain the absolute error in the scaled derivative structure factor h(kL) = 
(kLfg(kL): 



Ah{kL) = (kL) 2 Ag(kL) ~ 8 



1 



(26) 



_f tan" 1 (2/0 
This result indicates that the multi-step error in correlations increases with a 



5 Error in analytic time step 



In this Section, we calculate the error in analytic time step due to the error 
in scaled derivative structures (Ah) calculated in Section 4. 



5. 1 Conserved systems 



CH systems are purely dissipative systems — the energy density E monoton- 
ically decreases with the analytic time r with the relation E oc 1/L oc t -1 / 3 
[2]. We can calculate the analytic time in terms of the energy density E: 
t = Bch/E 3 , where the prefactor Bch can be numerically determined by 
requiring At = At as At — > in the late-time scaling regime since uncon- 
ditionally stable algorithms are arbitrarily accurate as At — > 0. We can then 
calculate the analytic time step by differentiating r with respect to E: 



Wc H AE -3A£t 4 / 3 



E 4 



B, 



1/3 
CH 



(27) 



and AE can be calculated by integrating from each Fourier mode: 



AE= I d 2 k 



ve 



6F 



(2tt) 2 U<ty fc , 



< )k (t + At) - <j> k (t)} 



d 2 k 



(2nk) 



-At eff (k,At)T k 



f28) 



where the time derivative dtfi-k/dr = —k 2 5F/S(pk from Eq. (J2J), and <pk{t + 
At) — (pkif) = At e ffd4>k/dr from Eq. ([6]) are used. The time-derivative corre- 
lation function Tf, [9"|27f28j has a natural scaling form given by 



T, = 



dr dr 




h(kL) 



L 2 h(kL) 
9r 4 /3 



(29) 



where L = Lqt 1 ^ 3 , h(x) is the 2D scaling function [28], and L is a constant. 
Taking the error in scaled structures into account, h(x) takes the form [28] 



h(x) 



a 



CH 



±Ah, 



x 



(30) 



as x > 1, where Cqh is a constant and Ah > is the absolute error in Eq. 
(12~TT) . The sign + or — will be determined later. Using the optimal driving 
scheme At = At 2 / 3 [23J, we can solve for AE in Eq. (|2"5|) and for the analytic 
time step At in Eq. (12 7\i . We have 



At: 



CcHL Q At 

,1/3 



VC 



ct.r 



<SxB%£ J x 2 [l + ( ai -l)Ax 2 /L 2 ] 



AhL 2 At 



m 



6ttB 



1/3 
CH 



dx 



x[l + (ai - l)Ax 2 /Ll] 



21 : 



(31) 



Solving the integral, we obtain 



Ar = At 



LoCchvOi 



d] vC4 + 0(A) 



(32) 



where 



a ~ — — In 



127TB. 



2/3 
CH 



T 2 



(tti - 1)4 



(33) 



and Eq. (12T|) is used. Note that we drop the + sign in Eq. (I3TI) according to 
the requirement of Ar < At implied by Eq. (J7|). We have 1 — At/ At ~ vA, 
which is consistent with the previous results [22J23J . 
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5.2 Non-conserved systems 



Similar to the calculation in conserved systems, we can calculate the analytic 
time step by differentiating r with respect to E: 



-2B AC AE -2Afir 3 / 2 

Ar = ~ P~ = -eg— (34) 

Integrating AE from each Fourier mode, we have 



1/5 1 K8F 



AE ~J d2k 7^{(wJ [Mt + At) - Mt) \ 



(2tt) 2 \\6<f> k/ 
1 



o 
1/5 

d 2 fc-^— At e/7 (/.-.An//,. !35i 



where the time derivative <90_fc/<9r = —5F/5(j)k from Eq. (j3j), and cj)k(t + At) — 
4>k(f) = Ateffd^k/dr from Eq. ([H]) are used. The time- derivative correlation 
function Tk [9.27.28] has a natural scaling form of 



where L = Lqt 1 / 2 , h(x) is the 2D scaling function [2B], and L is a constant. 
Taking the error in scaled structures into account, h(x) takes the form [28] 



h(x) = — ± Ah, (37) 

x 

as x > 1, where C^c is a constant and A/i > is the absolute error in Eq. 
(1261) . The sign + or — will be determined later. We can solve for AE in Eq. 
(1351) and for Ar from Eq. (J34J) . We have 



, i / / At(l-oa) 1 

C AC AtL tan W 1+At(ai-1) £ 

Ar = — 



,1/2 



4ttS^ J[l + At(a 1 -l)]At(l-a 2 ) 



AfeLgrV 2 / At(l - a 2 

:Sig(l -a 2 ) I + [1 + At( ai - 

As At = oo, we obtain the maximal analytic time step: 



± 8^(1 - o 2 ) ^ I 1 + [1 + A*(o! - l)]e) ' ( * 
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Ar 



ftan-Ha/O L 2 ln(l + aV^A/ir 1 / 2 



(ax - 1)(1 -o 2 ) 



87rBig(l 



(39) 



a 2 J 



where Eq. ( 1261) is used. Note that we drop the + sign in Eq. (1381) according to 
the requirement of Ar^ < l/(oi — 1) implied by Eq. 0. Define 



Atoo = Ar 00 A/(ai - 1)(1 - o 2 ) = £tan * 



f(a)r 



1/2 



(40) 



where 



/(a) 



T 2 



ttB 



1/2 
AC 



.£ tan" 1 (a/0 




(41) 




100 



Fig. 2. The discrepancy in the maximal analytic time step /(6)t ' 2 vs. the analytic 
time r in the scaling regime (20 < r < 100) with a = 2 (a\ = 3 and a 2 = —7) for 
system size L^ = 512. £ is taken to be 0.92 for best fit by eye. The discrepancy 
shows a t 1 / 2 behavior, as predicted in Eq. (|40p . 



In a previous paper [22], we fitted the time-averaged Ar^ with £tan _1 (a/£), 
and found that the fitting gets worse as a increases. This is consistent with 
the fact that /(a) monotonically increases with a. In addition, for a chosen a, 
in Fig. [21 we test the r-dependence of the discrepancy /(a)r 1//2 in the scaling 
regime and find excellent agreement. 
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6 Summary 



We find it is interesting to compare the analytic time r with the algorithmic 
time t. The analytic time r, determined by energy evolution, is a more funda- 
mental quantity. The unconditionally stable update Eq. (jBJ) is directly related 
to d<j)k/dr, and all scaling analysis is in analytic time r. The algorithmic time 
step At, on the other hand, is a more practical quantity. It appears in the driv- 
ing scheme and is used to express the Fourier space effective time step and to 
calculate the analytic time step. Our accuracy study on the scaled correlations 
g(x) and h(x) and subsequently on the analytic time step At based on the 
error accumulation relation Eq. (ITS]) were verified by numerical simulations. 
The analytic time step, the effective time step, the algorithmic time step, and 
the Fourier space multi step error are fundamental in understanding the er- 
ror behavior in unconditionally stable algorithms for coarsening systems. The 
methodology can be readily generalized to study other systems. We hope to 
report that work in a subsequent paper. 
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A Scaling of field derivatives in Fourier space 



In order to explore the accuracy of unconditionally stable algorithms in Fourier 
space, it is necessary to know the scaling of field derivatives for all modes. The 
structure factor S(k) = (\4>k\ 2 ) = L 2 g(kL), where g{kL) ~ (kL)~ 3 as k > 1/L 
[9]. Therefore we obtain 



<P k ~ L{kL)^' 2 = k-^L- 1 ' 2 . (A.l) 

Previous studies [91128] showed that dcfrk/dr = (dL/dr)k<pk as /cL > 1. Then 
we obtain the time-derivative correlation function 



T(k) = (\d<f) k /dT\ 2 ) = (dL/dT) 2 k 2 (\<j) k \ 2 ) = (dL/drfh^kL), (A.2) 

where the scaling function hi(kL) = k 2 L 2 g(kL) ~ (kL)" 1 ~ L~ l as k > 1/L. 
Subsequently we obtain 
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i fc dL 



^L^k- 1 ' 2 . (A.3) 

ar dr 

The generalization of higher order time-derivative correlations is (\d n 4>k/dr n \ 2 ) < 
(dL/dT) 2 k 2 (\d n ~ 1 4>k/dT n ~ 1 \ 2 ), where "~" indicates that generally the left hand 
side may not exactly equal to the right hand side. Applying this relation yields 

(\d n (f) k /dT n \ 2 ) ~ (dL/dT) 2n L 2 ~ 2n h n (kL), (A.4) 

where h n (kL) = k 2 L 2 h n _i(kL) ~ (kL) 2n ~ 3 as k > 1/L. Therefore we have 



dr n \dr J v ' 

-2n/3-i/6^,n-3/2 conserved 



n/2 i/4^,n 3/2 non-conserved 

The above expression is valid for n > for two dimensional scalar order 
parameter (s). 
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